# ===============================================
# Paris ratchet paper
# Code for extensions of randomization inference and sensitivity analysis
# ===============================================



#### Contents

#' 1. Load data for models (from "00_make_analysis_data.R")
#' 2. Randomization inference
#' 3. Sensitivity analysis


## Objects

#### Initialize #### 

library(tidyverse)
library(knitr)
library(DescTools) # Winsorize()
library(broom)
library(modelsummary)
library(magrittr) 
library(sensemakr) # for sensitivity 



# Set working directory


# Define core ISO3C codes to use in analysis
core_iso3c <- c("AFG", "AGO", "ALB", "AND", "ARE", "ARG", "ARM", "ATG", "AUS", "AUT", "AZE", "BDI", "BEL", "BEN", "BFA", "BGD", "BGR", "BHR", "BHS", "BIH", "BLR", "BLZ", "BOL", "BRA", "BRB", "BRN", "BTN", "BWA", "CAF", "CAN", "CHE", "CHL", "CHN", "CIV", "CMR", "COD", "COG", "COK", "COL", "COM", "CPV", "CRI", "CUB", "CYP", "CZE", "DEU", "DJI", "DMA", "DNK", "DOM", "DZA", "ECU", "EGY", "ERI", "ESP", "EST", "ETH", "FIN", "FJI", "FRA", "FSM", "GAB", "GBR", "GEO", "GHA", "GIN", "GMB", "GNB", "GNQ", "GRC", "GRD", "GTM", "GUY", "HND", "HRV", "HTI", "HUN", "IDN", "IND", "IRL", "IRN", "IRQ", "ISL", "ISR", "ITA", "JAM", "JOR", "JPN", "KAZ", "KEN", "KGZ", "KHM", "KIR", "KNA", "KOR", "KWT", "LAO", "LBN", "LBR", "LBY", "LCA", "LIE", "LKA", "LSO", "LTU", "LUX", "LVA", "MAR", "MDA", "MDG", "MDV", "MEX", "MHL", "MKD", "MLI", "MLT", "MMR", "MNE", "MNG", "MOZ", "MRT", "MUS", "MWI", "MYS", "NAM", "NER", "NGA", "NIC", "NLD", "NOR", "NPL", "NRU", "NZL", "OMN", "PAK", "PAN", "PER", "PHL", "PLW", "PNG", "POL", "PRK", "PRT", "PRY", "QAT", "ROU", "RUS", "RWA", "SAU", "SDN", "SEN", "SGP", "SLB", "SLE", "SLV", "SOM", "SRB", "SSD", "STP", "SUR", "SVK", "SVN", "SWE", "SWZ", "SYC", "SYR", "TCD", "TGO", "THA", "TJK", "TKM", "TLS", "TON", "TTO", "TUN", "TUR", "TUV", "TZA", "UGA", "UKR", "URY", "USA", "UZB", "VCT", "VEN", "VNM", "VUT", "WSM", "YEM", "ZAF", "ZMB", "ZWE")

data_for_models <- readRDS("output/data_for_models20230920.Rds")

#### Randomization inference ####

## Randomization inference
t1_a <- data_for_models %>% 
  lm(glasgow_target_safe ~ paris_target_safe + tradeflow_percentage, 
     data = .)
summary(t1_a)
t2_c <- data_for_models %>% 
  lm(glasgow_target_safe ~ paris_target_safe + io_percentage + tradeflow_percentage,
     data = .)
summary(t2_c)

data_for_ri <- data_for_models %>% 
  select(iso3c, paris_target_safe, glasgow_target_safe, tradeflow_percentage, io_percentage) %>% 
  filter(!is.na(paris_target_safe),
         !is.na(glasgow_target_safe))


set.seed(20230611)
seeds <- sample(1:1e6, size = 1000, replace = F)
store_beta <- rep(NA, times = length(seeds))
store_t <- rep(NA, times = length(seeds))

for (sim in seeds){
  # print(sim)
  set.seed(sim)
  draw_d <- sample(data_for_ri$io_percentage, size = 112, replace = T)
  model <- cbind.data.frame(draw_d, data_for_ri) %>%
    lm(glasgow_target_safe ~ paris_target_safe + tradeflow_percentage + 
         draw_d,
       data = .)
  store_beta[which(seeds==sim)] <- coefficients(model)["draw_d"]
  store_t[which(seeds==sim)] <- tidy(model) %>% filter(term == "draw_d") %>% pull(statistic)
}

as.data.frame(store_beta) %>% 
  select(beta = 1) %>% 
  reframe(beta, nobs = n()) %>% 
  # filter(abs(beta) > abs(coefficients(t2_c)[3])) %>% 
  filter(abs(beta) > abs(coefficients(t2_c)[3])) %>% 
  reframe(beta, nobs, remaining = n(), ri_p = remaining/nobs) 
# Only 2-3 times per 1000 is the placebo treatment larger than the actual treatment
# In 100,000 trials, only 195 larger than ATE (0.00195) (.195%)

ri_plot <- as.data.frame(store_beta) %>% 
  ggplot(., aes(store_beta)) +
  geom_histogram(aes(y = ..density..),
                 bins = 50,
                 fill = "grey20",
                 color = "grey90") +
  geom_vline(xintercept = coefficients(t2_c)["io_percentage"], linetype = 2, color = "red") +
  geom_vline(xintercept = -coefficients(t2_c)["io_percentage"], linetype = 3, color = "red") +
  labs(color = "", shape = "", 
       x = "Estimated ATE", y = "Density") +
  theme(legend.position = "bottom", 
        panel.background = element_rect(fill=NA),
        panel.grid.major.x = element_blank(),
        panel.grid.major.y = element_blank(),
        panel.grid.minor.y =  element_blank(),
        panel.grid.minor.x =  element_blank(),
        axis.line = element_line(size = 1, colour = "black"),
        axis.text = element_text(size = 12),
        strip.text = element_text(size = 12),
        legend.text = element_text(size = 12),
        legend.title = element_text(size = 12),
        axis.title = element_text(size = 12),
        panel.ontop = FALSE)
ggsave(plot = ri_plot, 
       filename = "plots/ri_W_IO.pdf", units = 'in', height = 6, width = 8)






#### Sensitivity analysis #### 

sen1 <- data_for_models %>% 
  rename(`Paris target` = paris_target_safe) %>% 
  lm(glasgow_target_safe ~ `Paris target` + io_percentage + tradeflow_percentage,
     data = .)


sensitivity_paris <- sensemakr(model = sen1, 
                               treatment = "io_percentage",
                               benchmark_covariates = "`Paris target`",
                               kd = 1)
print(sensitivity_paris)
plot(sensitivity_paris)
summary(sensitivity_paris)
ovb_minimal_reporting(sensitivity_paris)

#' Partial R2 of treatment with outcome: 0.0426 
#' Robustness value is 0.1898
#' Bound on R2dz.x = 0.0122
#' Bound on R2yz.dx = 0.71 
#' In plot(), a confounder as strong as paris_target_safe, 
#' would reduce value of tradeflow_percentage by 55%
#' From 3.3 to 1.814
#' Bound of R2dz.x is lower than 'Partial R2 of treatment with outcome'
#' so even extreme confounder explaining all residual variation of outcome
#' and as strong associated with treatment as paris_target_safe would not reduce coefficient to zero

pdf(file = "plots/sens_plot.pdf")
plot(sensitivity_paris)
dev.off()
